#!/usr/bin/env python
"""\
NAME
    mcubes.py - Marching cubes contouring of 3D data

CUBE REFERENCE:
   Vertices                   Edges
       4             7                  8         
       +-------------+           +-------------+  
      /|            /|        5 /|           7/|  
    5/ |         6 / |         / |    6      / |  
    +-------------+  |        +-------------+  |  
    |  |          |  |        | 9|          |  |12 
    |  |0         | 3|     10 |  |    4   11|  |  
    |  +----------|--+        |  +----------|--+  
    | /           | /         | /1          | /3  
    |/1          2|/          |/     2      |/    
    +-------------+           +-------------+     

From p. 216 of Watt and Watt.

Copyright (c) 2003 Richard P. Muller (rmuller@sandia.gov). All rights
reserved. See the LICENSE file for licensing details.
"""

import os,profile,pstats

# Returns the vertices that make up an edge
edge2pts = [(0,1),(1,2),(2,3),(3,0),(4,5),(5,6),(6,7),(7,4),
            (0,4),(1,5),(2,6),(3,7)]

# For a given case (0-255), this array yields the corresponding
# triangles, given as a tuple of three edges
tris_as_edges = [
    [],
    [(0, 8, 3)],
    [(0, 1, 9)],
    [(1, 8, 3), (9, 8, 1)],
    [(1, 2, 10)],
    [(0, 8, 3), (1, 2, 10)],
    [(10, 9, 0), (0, 2, 10)],
    [(2, 8, 3), (2, 10, 8), (10, 9, 8)],
    [(3, 11, 2)],
    [(2, 0, 8), (8, 11, 2)],
    [(1, 9, 0), (2, 3, 11)],
    [(1, 11, 2), (1, 9, 11), (9, 8, 11)],
    [(3, 11, 1), (1, 11, 10)],
    [(0, 10, 1), (0, 8, 10), (11, 10, 8)],
    [(3, 9, 0), (3, 11, 9), (11, 10, 9)],
    [(9, 8, 11), (10, 9, 11)],
    [(4, 7, 8)],
    [(4, 3, 0), (7, 3, 4)],
    [(0, 1, 9), (8, 4, 7)],
    [(4, 1, 9), (4, 7, 1), (7, 3, 1)],
    [(1, 2, 10), (8, 4, 7)],
    [(3, 4, 7), (3, 0, 4), (1, 2, 10)],
    [(9, 2, 10), (9, 0, 2), (8, 4, 7)],
    [(2, 10, 9), (2, 9, 7), (2, 7, 3), (7, 9, 4)],
    [(8, 4, 7), (3, 11, 2)],
    [(11, 4, 7), (11, 2, 4), (2, 0, 4)],
    [(9, 0, 1), (8, 4, 7), (2, 3, 11)],
    [(4, 7, 11), (9, 4, 11), (9, 11, 2), (2, 1, 9)],
    [(3, 11, 10), (10, 1, 3), (7, 8, 4)],
    [(11, 10, 1), (4, 11, 1), (7, 11, 4), (4, 1, 0)],
    [(4, 7, 8), (3, 11, 9), (9, 11, 10), (9, 0, 3)],
    [(4, 7, 11), (4, 11, 9), (9, 11, 10)],
    [(9, 5, 4)],
    [(9, 5, 4), (0, 8, 3)],
    [(0, 5, 4), (1, 5, 0)],
    [(8, 5, 4), (8, 3, 5), (3, 1, 5)],
    [(1, 2, 10), (9, 5, 4)],
    [(3, 0, 8), (1, 2, 10), (4, 9, 5)],
    [(5, 2, 10), (5, 4, 2), (4, 0, 2)],
    [(2, 10, 5), (3, 2, 5), (3, 5, 4), (3, 4, 8)],
    [(9, 5, 4), (2, 3, 11)],
    [(0, 11, 2), (0, 8, 11), (4, 9, 5)],
    [(0, 5, 4), (0, 1, 5), (2, 3, 11)],
    [(2, 1, 5), (2, 5, 8), (2, 8, 11), (4, 8, 5)],
    [(10, 3, 11), (10, 1, 3), (9, 5, 4)],
    [(4, 9, 5), (0, 8, 1), (8, 10, 1), (8, 11, 10)],
    [(5, 4, 0), (5, 0, 11), (5, 11, 10), (11, 0, 3)],
    [(5, 4, 8), (5, 8, 10), (10, 8, 11)],
    [(9, 7, 8), (5, 7, 9)],
    [(9, 3, 0), (9, 5, 3), (5, 7, 3)],
    [(0, 7, 8), (0, 1, 7), (1, 5, 7)],
    [(1, 5, 3), (3, 5, 7)],
    [(9, 7, 8), (9, 5, 7), (10, 1, 2)],
    [(10, 1, 2), (9, 5, 0), (5, 3, 0), (5, 7, 3)],
    [(8, 0, 2), (8, 2, 5), (8, 5, 7), (10, 5, 2)],
    [(2, 10, 5), (2, 5, 3), (3, 5, 7)],
    [(7, 9, 5), (7, 8, 9), (3, 11, 2)],
    [(9, 5, 7), (9, 7, 2), (9, 2, 0), (2, 7, 11)],
    [(2, 3, 11), (0, 1, 8), (1, 7, 8), (1, 5, 7)],
    [(11, 2, 1), (11, 1, 7), (7, 1, 5)],
    [(9, 5, 8), (8, 5, 7), (10, 1, 3), (10, 3, 11)],
    [(5, 7, 0), (5, 0, 9), (7, 11, 0), (1, 0, 10), (11, 10, 0)],
    [(11, 10, 0), (11, 0, 3), (10, 5, 0), (8, 0, 7), (5, 7, 0)],
    [(5, 11, 10), (11, 5, 7)],
    [(10, 6, 5)],
    [(0, 8, 3), (6, 5, 10)],
    [(9, 0, 1), (5, 10, 6)],
    [(1, 8, 3), (1, 9, 8), (5, 10, 6)],
    [(1, 6, 5), (2, 6, 1)],
    [(1, 6, 5), (1, 2, 6), (3, 0, 8)],
    [(9, 6, 5), (9, 0, 6), (0, 2, 6)],
    [(5, 9, 8), (5, 8, 2), (5, 2, 6), (3, 2, 8)],
    [(2, 3, 11), (10, 6, 5)],
    [(11, 0, 8), (11, 2, 0), (10, 6, 5)],
    [(0, 1, 9), (2, 3, 11), (5, 10, 6)],
    [(5, 10, 6), (1, 9, 2), (9, 11, 2), (9, 8, 11)],
    [(6, 3, 11), (6, 5, 3), (5, 1, 3)],
    [(0, 8, 11), (0, 11, 5), (0, 5, 1), (5, 11, 6)],
    [(3, 11, 6), (0, 3, 6), (0, 6, 5), (0, 5, 9)],
    [(6, 5, 9), (6, 9, 11), (11, 9, 8)],
    [(5, 10, 6), (4, 7, 8)],
    [(4, 3, 0), (4, 7, 3), (6, 5, 10)],
    [(1, 9, 0), (5, 10, 6), (8, 4, 7)],
    [(10, 6, 5), (1, 9, 7), (1, 7, 3), (7, 9, 4)],
    [(6, 1, 2), (6, 5, 1), (4, 7, 8)],
    [(1, 2, 5), (5, 2, 6), (3, 0, 4), (3, 4, 7)],
    [(8, 4, 7), (9, 0, 5), (0, 6, 5), (0, 2, 6)],
    [(7, 3, 9), (7, 9, 4), (3, 2, 9), (5, 9, 6), (2, 6, 9)],
    [(3, 11, 2), (7, 8, 4), (10, 6, 5)],
    [(5, 10, 6), (4, 7, 2), (4, 2, 0), (2, 7, 11)],
    [(0, 1, 9), (4, 7, 8), (2, 3, 11), (5, 10, 6)],
    [(9, 2, 1), (9, 11, 2), (9, 4, 11), (7, 11, 4), (5, 10, 6)],
    [(8, 4, 7), (3, 11, 5), (3, 5, 1), (5, 11, 6)],
    [(5, 1, 11), (5, 11, 6), (1, 0, 11), (7, 11, 4), (0, 4, 11)],
    [(0, 5, 9), (0, 6, 5), (0, 3, 6), (11, 6, 3), (8, 4, 7)],
    [(6, 5, 9), (6, 9, 11), (4, 7, 9), (7, 11, 9)],
    [(6, 4, 9), (9, 10, 6)],
    [(4, 10, 6), (4, 9, 10), (0, 8, 3)],
    [(10, 0, 1), (10, 6, 0), (6, 4, 0)],
    [(8, 3, 1), (8, 1, 6), (8, 6, 4), (6, 1, 10)],
    [(1, 4, 9), (1, 2, 4), (2, 6, 4)],
    [(3, 0, 8), (1, 2, 9), (2, 4, 9), (2, 6, 4)],
    [(0, 2, 4), (4, 2, 6)],
    [(8, 3, 2), (8, 2, 4), (4, 2, 6)],
    [(10, 4, 9), (10, 6, 4), (11, 2, 3)],
    [(0, 8, 2), (2, 8, 11), (4, 9, 10), (4, 10, 6)],
    [(3, 11, 2), (0, 1, 6), (0, 6, 4), (6, 1, 10)],
    [(6, 4, 1), (6, 1, 10), (4, 8, 1), (2, 1, 11), (8, 11, 1)],
    [(9, 6, 4), (9, 3, 6), (9, 1, 3), (11, 6, 3)],
    [(8, 11, 1), (8, 1, 0), (11, 6, 1), (9, 1, 4), (6, 4, 1)],
    [(3, 11, 6), (3, 6, 0), (0, 6, 4)],
    [(6, 4, 8), (8, 11, 6)],
    [(7, 10, 6), (7, 8, 10), (8, 9, 10)],
    [(0, 7, 3), (0, 10, 7), (0, 9, 10), (6, 7, 10)],
    [(10, 6, 7), (1, 10, 7), (1, 7, 8), (1, 8, 0)],
    [(10, 6, 7), (10, 7, 1), (1, 7, 3)],
    [(1, 2, 6), (1, 6, 8), (1, 8, 9), (8, 6, 7)],
    [(2, 6, 9), (2, 9, 1), (6, 7, 9), (0, 9, 3), (7, 3, 9)],
    [(7, 8, 0), (7, 0, 6), (6, 0, 2)],
    [(7, 3, 2), (6, 7, 2)],
    [(2, 3, 11), (10, 6, 8), (10, 8, 9), (8, 6, 7)],
    [(2, 0, 7), (2, 7, 11), (0, 9, 7), (6, 7, 10), (9, 10, 7)],
    [(1, 8, 0), (1, 7, 8), (1, 10, 7), (6, 7, 10), (2, 3, 11)],
    [(11, 2, 1), (11, 1, 7), (10, 6, 1), (6, 7, 1)],
    [(8, 9, 6), (8, 6, 7), (9, 1, 6), (11, 6, 3), (1, 3, 6)],
    [(0, 9, 1), (11, 6, 7)],
    [(7, 8, 0), (7, 0, 6), (3, 11, 0), (11, 6, 0)],
    [(7, 11, 6)],
    [(7, 6, 11)],
    [(3, 0, 8), (11, 7, 6)],
    [(0, 1, 9), (11, 7, 6)],
    [(8, 1, 9), (8, 3, 1), (11, 7, 6)],
    [(10, 1, 2), (6, 11, 7)],
    [(1, 2, 10), (3, 0, 8), (6, 11, 7)],
    [(2, 9, 0), (2, 10, 9), (6, 11, 7)],
    [(6, 11, 7), (2, 10, 3), (10, 8, 3), (10, 9, 8)],
    [(7, 2, 3), (6, 2, 7)],
    [(7, 0, 8), (7, 6, 0), (6, 2, 0)],
    [(2, 7, 6), (2, 3, 7), (0, 1, 9)],
    [(1, 6, 2), (1, 8, 6), (1, 9, 8), (8, 7, 6)],
    [(10, 7, 6), (10, 1, 7), (1, 3, 7)],
    [(10, 7, 6), (1, 7, 10), (1, 8, 7), (1, 0, 8)],
    [(0, 3, 7), (0, 7, 10), (0, 10, 9), (6, 10, 7)],
    [(7, 6, 10), (7, 10, 8), (8, 10, 9)],
    [(8, 4, 6), (6, 11, 8)],
    [(3, 6, 11), (3, 0, 6), (0, 4, 6)],
    [(8, 6, 11), (8, 4, 6), (9, 0, 1)],
    [(9, 4, 6), (9, 6, 3), (9, 3, 1), (11, 3, 6)],
    [(6, 8, 4), (6, 11, 8), (2, 10, 1)],
    [(1, 2, 10), (3, 0, 11), (0, 6, 11), (0, 4, 6)],
    [(4, 11, 8), (4, 6, 11), (0, 2, 9), (2, 10, 9)],
    [(10, 9, 3), (10, 3, 2), (9, 4, 3), (11, 3, 6), (4, 6, 3)],
    [(8, 2, 3), (8, 4, 2), (4, 6, 2)],
    [(0, 4, 2), (4, 6, 2)],
    [(1, 9, 0), (2, 3, 4), (2, 4, 6), (4, 3, 8)],
    [(1, 9, 4), (1, 4, 2), (2, 4, 6)],
    [(8, 1, 3), (8, 6, 1), (8, 4, 6), (6, 10, 1)],
    [(10, 1, 0), (10, 0, 6), (6, 0, 4)],
    [(4, 6, 3), (4, 3, 8), (6, 10, 3), (0, 3, 9), (10, 9, 3)],
    [(10, 9, 4), (4, 6, 10)],
    [(4, 9, 5), (7, 6, 11)],
    [(0, 8, 3), (4, 9, 5), (11, 7, 6)],
    [(5, 0, 1), (5, 4, 0), (7, 6, 11)],
    [(11, 7, 6), (8, 3, 4), (3, 5, 4), (3, 1, 5)],
    [(9, 5, 4), (10, 1, 2), (7, 6, 11)],
    [(6, 11, 7), (1, 2, 10), (0, 8, 3), (4, 9, 5)],
    [(7, 6, 11), (5, 4, 10), (4, 2, 10), (4, 0, 2)],
    [(3, 4, 8), (3, 5, 4), (3, 2, 5), (10, 5, 2), (11, 7, 6)],
    [(7, 2, 3), (7, 6, 2), (5, 4, 9)],
    [(9, 5, 4), (0, 8, 6), (0, 6, 2), (6, 8, 7)],
    [(3, 6, 2), (3, 7, 6), (1, 5, 0), (5, 4, 0)],
    [(6, 2, 8), (6, 8, 7), (2, 1, 8), (4, 8, 5), (1, 5, 8)],
    [(9, 5, 4), (10, 1, 6), (1, 7, 6), (1, 3, 7)],
    [(1, 6, 10), (1, 7, 6), (1, 0, 7), (8, 7, 0), (9, 5, 4)],
    [(4, 0, 10), (4, 10, 5), (0, 3, 10), (6, 10, 7), (3, 7, 10)],
    [(7, 6, 10), (7, 10, 8), (5, 4, 10), (4, 8, 10)],
    [(6, 9, 5), (6, 11, 9), (11, 8, 9)],
    [(3, 6, 11), (0, 6, 3), (0, 5, 6), (0, 9, 5)],
    [(0, 11, 8), (0, 5, 11), (0, 1, 5), (5, 6, 11)],
    [(6, 11, 3), (6, 3, 5), (5, 3, 1)],
    [(1, 2, 10), (9, 5, 11), (9, 11, 8), (11, 5, 6)],
    [(0, 11, 3), (0, 6, 11), (0, 9, 6), (5, 6, 9), (1, 2, 10)],
    [(11, 8, 5), (11, 5, 6), (8, 0, 5), (10, 5, 2), (0, 2, 5)],
    [(6, 11, 3), (6, 3, 5), (2, 10, 3), (10, 5, 3)],
    [(5, 8, 9), (5, 2, 8), (5, 6, 2), (3, 8, 2)],
    [(9, 5, 6), (9, 6, 0), (0, 6, 2)],
    [(1, 5, 8), (1, 8, 0), (5, 6, 8), (3, 8, 2), (6, 2, 8)],
    [(1, 5, 6), (2, 1, 6)],
    [(1, 3, 6), (1, 6, 10), (3, 8, 6), (5, 6, 9), (8, 9, 6)],
    [(10, 1, 0), (10, 0, 6), (9, 5, 0), (5, 6, 0)],
    [(0, 3, 8), (5, 6, 10)],
    [(5, 6, 10)], 
    [(11, 5, 10), (7, 5, 11)],
    [(11, 5, 10), (11, 7, 5), (8, 3, 0)],
    [(5, 11, 7), (5, 10, 11), (1, 9, 0)],
    [(10, 7, 5), (10, 11, 7), (9, 8, 1), (8, 3, 1)],
    [(11, 1, 2), (11, 7, 1), (7, 5, 1)],
    [(0, 8, 3), (1, 2, 7), (1, 7, 5), (7, 2, 11)],
    [(9, 7, 5), (9, 2, 7), (9, 0, 2), (2, 11, 7)],
    [(7, 5, 2), (7, 2, 11), (5, 9, 2), (3, 2, 8), (9, 8, 2)],
    [(2, 5, 10), (2, 3, 5), (3, 7, 5)],
    [(8, 2, 0), (8, 5, 2), (8, 7, 5), (10, 2, 5)],
    [(9, 0, 1), (5, 10, 3), (5, 3, 7), (3, 10, 2)],
    [(9, 8, 2), (9, 2, 1), (8, 7, 2), (10, 2, 5), (7, 5, 2)],
    [(1, 3, 5), (3, 7, 5)],
    [(0, 8, 7), (0, 7, 1), (1, 7, 5)],
    [(9, 0, 3), (9, 3, 5), (5, 3, 7)],
    [(9, 8, 7), (5, 9, 7)],
    [(5, 8, 4), (5, 10, 8), (10, 11, 8)],
    [(5, 0, 4), (5, 11, 0), (5, 10, 11), (11, 3, 0)],
    [(0, 1, 9), (8, 4, 10), (8, 10, 11), (10, 4, 5)],
    [(10, 11, 4), (10, 4, 5), (11, 3, 4), (9, 4, 1), (3, 1, 4)],
    [(2, 5, 1), (2, 8, 5), (2, 11, 8), (4, 5, 8)],
    [(0, 4, 11), (0, 11, 3), (4, 5, 11), (2, 11, 1), (5, 1, 11)],
    [(0, 2, 5), (0, 5, 9), (2, 11, 5), (4, 5, 8), (11, 8, 5)],
    [(9, 4, 5), (2, 11, 3)],
    [(2, 5, 10), (3, 5, 2), (3, 4, 5), (3, 8, 4)],
    [(5, 10, 2), (5, 2, 4), (4, 2, 0)],
    [(3, 10, 2), (3, 5, 10), (3, 8, 5), (4, 5, 8), (0, 1, 9)],
    [(5, 10, 2), (5, 2, 4), (1, 9, 2), (9, 4, 2)],
    [(8, 4, 5), (8, 5, 3), (3, 5, 1)],
    [(0, 4, 5), (1, 0, 5)],
    [(8, 4, 5), (8, 5, 3), (9, 0, 5), (0, 3, 5)],
    [(9, 4, 5)],
    [(4, 11, 7), (4, 9, 11), (9, 10, 11)],
    [(0, 8, 3), (4, 9, 7), (9, 11, 7), (9, 10, 11)],
    [(1, 10, 11), (1, 11, 4), (1, 4, 0), (7, 4, 11)],
    [(3, 1, 4), (3, 4, 8), (1, 10, 4), (7, 4, 11), (10, 11, 4)],
    [(4, 11, 7), (9, 11, 4), (9, 2, 11), (9, 1, 2)],
    [(9, 7, 4), (9, 11, 7), (9, 1, 11), (2, 11, 1), (0, 8, 3)],
    [(11, 7, 4), (11, 4, 2), (2, 4, 0)],
    [(11, 7, 4), (11, 4, 2), (8, 3, 4), (3, 2, 4)],
    [(2, 9, 10), (2, 7, 9), (2, 3, 7), (7, 4, 9)],
    [(9, 10, 7), (9, 7, 4), (10, 2, 7), (8, 7, 0), (2, 0, 7)],
    [(3, 7, 10), (3, 10, 2), (7, 4, 10), (1, 10, 0), (4, 0, 10)],
    [(1, 10, 2), (8, 7, 4)],
    [(4, 9, 1), (4, 1, 7), (7, 1, 3)],
    [(4, 9, 1), (4, 1, 7), (0, 8, 1), (8, 7, 1)],
    [(4, 0, 3), (7, 4, 3)],
    [(4, 8, 7)],
    [(9, 11, 8), (11, 9, 10)],
    [(3, 0, 9), (3, 9, 11), (11, 9, 10)],
    [(0, 1, 10), (0, 10, 8), (8, 10, 11)],
    [(3, 1, 11), (11, 1, 10)],
    [(1, 2, 11), (1, 11, 9), (9, 11, 8)],
    [(3, 0, 9), (3, 9, 11), (1, 2, 9), (2, 11, 9)],
    [(0, 2, 11), (11, 8, 0)],
    [(3, 2, 11)],
    [(2, 3, 8), (2, 8, 10), (10, 8, 9)],
    [(2, 0, 9), (9, 10, 2)],
    [(2, 3, 8), (2, 8, 10), (0, 1, 8), (1, 10, 8)],
    [(10, 2, 1)],
    [(1, 3, 8), (9, 1, 8)],
    [(0, 9, 1)],
    [(3, 8, 0)],
    [],
    ]

def which_case((v0,v1,v2,v3,v4,v5,v6,v7),contour_value):
    # should rewrite using bitarrays
    value = 0
    if v0 > contour_value: value += 1
    if v1 > contour_value: value += 2
    if v2 > contour_value: value += 4
    if v3 > contour_value: value += 8
    if v4 > contour_value: value += 16
    if v5 > contour_value: value += 32
    if v6 > contour_value: value += 64
    if v7 > contour_value: value += 128
    return value

def contour(contour_value,
            p1,p2,p3,p4,p5,p6,p7,p8,
            g1,g2,g3,g4,g5,g6,g7,g8):
    "Return a list of triangles for this cube"
    # p1-p8 are tuples of xyz values for the vertices
    # g1-g8 are the corresponding magnitudes on the grid 
    tris = []
    case = which_case((g1,g2,g3,g4,g5,g6,g7,g8),contour_value)
    for e1,e2,e3 in tris_as_edges[case]:
        #xyz1 = edge_coordinates(e1,(p1,p2,p3,p4,p5,p6,p7,p8))
        #xyz2 = edge_coordinates(e2,(p1,p2,p3,p4,p5,p6,p7,p8))
        #xyz3 = edge_coordinates(e3,(p1,p2,p3,p4,p5,p6,p7,p8))
        xyz1 = interp_edge_coordinates(e1,(p1,p2,p3,p4,p5,p6,p7,p8),
                                       (g1,g2,g3,g4,g5,g6,g7,g8),
                                       contour_value)
        xyz2 = interp_edge_coordinates(e2,(p1,p2,p3,p4,p5,p6,p7,p8),
                                       (g1,g2,g3,g4,g5,g6,g7,g8),
                                       contour_value)
        xyz3 = interp_edge_coordinates(e3,(p1,p2,p3,p4,p5,p6,p7,p8),
                                       (g1,g2,g3,g4,g5,g6,g7,g8),
                                       contour_value)
        tris.append((xyz1,xyz2,xyz3))
    return tris

def edge_coordinates(edge,pts):
    start,end = edge2pts[edge]
    x0,y0,z0 = pts[start]
    x1,y1,z1 = pts[end]
    return 0.5*(x0+x1),0.5*(y0+y1),0.5*(z0+z1)

def interp_edge_coordinates(edge,pts,vals,contour_value):
    start,end = edge2pts[edge]
    vstart = vals[start]
    vend = vals[end]
    frac = (contour_value-vstart)/(vend-vstart)
    cfrac = 1-frac #conjugate of frac
    x0,y0,z0 = pts[start]
    x1,y1,z1 = pts[end]
    return cfrac*x0+frac*x1,cfrac*y0+frac*y1,cfrac*z0+frac*z1


def test_contour_sphere(npts):
    "Test example using POVRay module"
    from numpy import arange,array
    from Pistol.POVRay import Scene,TriangleMesh

    def sphere(x,y,z): return x*x+y*y+z*z

    # The goal here was to see whether MC would work when there weren't cubes,
    # e.g. when yvec = 0.1,1.0,0.0. It doesn't appear to.

    origin = array([-5.0,-5.0,-5.0])
    length = 10
    xvec = length*array([1.0,0.0,0.0])
    yvec = length*array([0.0,1.0,0.0])
    zvec = length*array([0.0,0.0,1.0])

    nx=ny=nz=npts

    triangles = []
    for i in range(nx-1):
        ifrac = i/(nx-1.)
        ifracp = (i+1.)/(nx-1.)
        for j in range(ny-1):
            jfrac = j/(ny-1.)
            jfracp = (j+1.)/(ny-1.)
            for k in range(nz-1):
                kfrac = k/(nz-1.)
                kfracp = (k+1.)/(nz-1.)
                x,y,z = origin + ifrac*xvec + jfrac*yvec + kfrac*zvec
                xp,yp,zp = origin + ifracp*xvec + jfracp*yvec + kfracp*zvec

                newtris = contour(3.0,
                                  (x,y,z),(x,y,zp),(x,yp,zp),(x,yp,z),
                                  (xp,y,z),(xp,y,zp),(xp,yp,zp),(xp,yp,z),
                                  sphere(x,y,z),sphere(x,y,zp),
                                  sphere(x,yp,zp),sphere(x,yp,z),
                                  sphere(xp,y,z),sphere(xp,y,zp),
                                  sphere(xp,yp,zp),sphere(xp,yp,z))
                if newtris:
                    for tri in newtris: triangles.append(tri)
    pov = Scene('test_sphere')
    tmesh = TriangleMesh(triangles,(0,0,1))
    pov.add(tmesh)
    pov.write_pov()
    pov.render()
    pov.display()
    return




